#Code built for Gabriel Fulton thesis in order to import raw text files
#Imported from text for processing, graphing, and analysis
#packages: quantmod, DT (for seperate tables)
#packages?: TTR, xts, zoo

#Set wd to whichever folder currently has the files which will be analyzed
#Surface
setwd("C:/Users/gfult/Dropbox/Salt/InProgress/Spectro Analysis/08Jun2018/BeetSpectraRedux/PSR+3500_1596061/2018_Jun_08")
#Home PC
#setwd("C:/Users/Gabriel Fulton/Dropbox/Salt/InProgress/Spectro Analysis/08Jun2018/BeetSpectraRedux/PSR+3500_1596061/2018_Jun_08")

#Read from .sed file based on spacing using fwf, "View(anysample)" for check
S1=read.fwf("Jun8beet_00031.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S2=read.fwf("Jun8beet_00032.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S3=read.fwf("Jun8beet_00033.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S4=read.fwf("Jun8beet_00034.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S5=read.fwf("Jun8beet_00035.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S6=read.fwf("Jun8beet_00036.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S7=read.fwf("Jun8beet_00037.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S8=read.fwf("Jun8beet_00038.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S9=read.fwf("Jun8beet_00039.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S10=read.fwf("Jun8beet_00040.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)

#Graph all individual readings (Currently can't fit 10...) [layout function cool af]
#x11() makes a serpate window for graph
#x11()
#par(mfrow=c(2,5))
#plot(S1)
#plot(S2)
#plot(S3)
#plot(S4)
#plot(S5)
#plot(S6)
#plot(S7)
#plot(S8)
#plot(S9)
#plot(S10)


#Average these spectral readings with each other
beetdish10=(S1+S2+S3+S4+S5+S6+S7+S8+S9+S10)/10

#Plot the cumulative average
x11()
par(mfrow=c(1,1))
plot(beetdish10, main="")

#Divided by highest value such that max is now 1
reflectance=beetdish10$Reflectance/max(beetdish10$Reflectance)
#x11()
#par(mfrow=c(1,1))

beetdish10DividedByMax=cbind.data.frame(350:2500,reflectance)
#plot(beetdishsalt5DividedByMax)

#Correlation
#Truncating data for correlation, 350-1400
#Analysis of key points (Salt) - 1490, 1990
#this covers 350 to 1400 nm wavelength
tbeetdish10=beetdish10[1:1051,1:2]
cor(tsalt$Reflectance,tbeetdish10$Reflectance)

#Divide out the water and dish
cor(tsalt$Reflectance,tbeetdish10$Reflectance/twaterdish$Reflectance)

#Taking natural log
#maybe divide out beet juice equivalence instead of waterdish0
cor(log(tsalt$Reflectance),log(tbeetdish10$Reflectance/twaterdish$Reflectance))

#Taking natural log as step 1
#maybe divide out beet juice equivalence instead of waterdish0
cor(log(tsalt$Reflectance),log(tbeetdish10$Reflectance)/log(twaterdish$Reflectance))

#Grab desired POI y values (at wavelength = 986,1001,1030,1058,1193,1256,1457,1632)
POIbeetdish10=rbind(beetdish10[637,],beetdish10[652,],beetdish10[681,],beetdish10[709,],beetdish10[844,],beetdish10[907,],beetdish10[1108,],beetdish10[1283,])
POIbeetdish10

#Salt derived POI
POI2beetdish10=rbind(beetdish10[129,],beetdish10[509,],beetdish10[557,],beetdish10[637,],beetdish10[744,],beetdish10[874,],beetdish10[921,],beetdish10[1308,])
POI2beetdish10

#Prints to notepad to copy to excel
write.table(beetdish10,"C:/Users/gfult/Desktop/Upload/beetdish10.txt",sep="\t")

#EVERYTHING BELOW HERE IS WORK IN PROGRESS FOR AUTOMATION OF PEAK FINDING

#Apparently embedded findPeaks is no good so adding this function:
find_peaks <- function (x, m = 3){
  shape <- diff(sign(diff(x, na.pad = FALSE)))
  pks <- sapply(which(shape < 0), FUN = function(i){
    z <- i - m + 1
    z <- ifelse(z > 0, z, 1)
    w <- i + m + 1
    w <- ifelse(w < length(x), w, length(x))
    if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
  })
  pks <- unlist(pks)
  pks
}
#Where m is number of points on either side of the peak
#the function can also be used to find local minima of any sequential vector x via find_peaks(-x)

#THIS AREA IS MESSED UP AND NEEDS WORK
#Peak and valley analysis (NOT SURE OF WHAT m IS APPROPRIATE, originally 3)
peaks=find_peaks(Stotal$Reflectance, m=30)
valleys=find_peaks(-Stotal$Reflectance, m=30)

#peaks=findPeaks(Stotal$Reflectance,thresh=0.05)
#valleys=findValleys(Stotal$Reflectance,thresh=0.05)

#Finds corresponding reflectance values for found wavelength peaks and valleys
ppeaks=Stotal[Stotal$Wavelength%in%peaks,]
pvalleys=Stotal[Stotal$Wavelength%in%valleys,]


#Highlight of said peak and valleys with corresponding wavelengths/severity
points(ppeaks,col="red",pch=19)
points(pvalleys,col="blue",pch=19)

#Creation of tables of corresponding data:
#Peaks:
View(ppeaks, "Peaks")
#Valleys:
View(pvalleys, "Valleys")

#For Table writing:
#write.xlsx(x, file, sheetName="Sheet1")

